In emission tomography the density data of particle emissions (in the following: emission density data) in a volume has to be found on the basis of measurement data. Emitting particles are for instance positrons or gamma photons. Emitted particles may get scattered or get absorbed, and may also change their type in the volume or before detection. Emitted particles or the particles originating from the emitted particles arrive in detector devices that report detection events. The above mentioned physical phenomena complicate the tomographic reconstruction, and for an appropriate result these have to be taken into account in the course of the reconstruction. An emission tomographic reconstruction method reconstructs the emission density data based on the measurement data obtained from signals of detector devices comprising a plurality of detector elements.
In Positron Emission Tomography (PET) the density of radioactive tracer materials has to be found, which is proportional with its emission density data. In the case of PET, the particles emitted from the tracer material are positrons, which are replaced by pairs of approximately oppositely directed gamma photons when they annihilate with electrons (cf. Physics reference manual, Geant4 9.1. Technical Report, CERN, 2007). These gamma photons may get scattered in the body and/or in the detector elements before they get finally absorbed in a detector element comprised in the detector device or module. The detector elements are usually in the form of detector crystals. In PET the detector device should detect pairs of gamma photons arriving almost simultaneously, so an event occurs if two detector elements nearly simultaneously detect two photons. This means that in PET, a pair of detector elements, sometimes called detector crystals, represents a conceptual detector, which is also called Line Of Response (LOR).
In Single Photon Emission Tomography (SPECT), emitted particles are gamma photons. A particle may also get scattered or absorbed in the object, in the collimator or in the detector crystals comprised in a detector device rotating around the object to be investigated. A detector crystal—i.e. a detector element—is able to detect gamma photons that went through a collimator associated to it. The collimator can be for example parallel hole, pinhole, multi-pinhole, fan or cone beam, slit-slat type or any other type. Thus in SPECT, a conceptual detector is defined by the actual orientation of the detector device, i.e. of its detector elements and the corresponding collimators. The line characterizing this orientation is called the projection line.
Two main types of the reconstruction methods are known. If the emission density is reconstructed only from the number of events per conceptual detector, then the reconstruction process is said to be based on binned data. In a binned type reconstruction the inputs comprise the total number of hits in conceptual detectors acquired during the measurement. Consequently, when binned reconstruction is applied, individual detector events undergo binning and the emission density data has to be reconstructed from the spatial histogram of photon pair hits in the case of PET or from the photon hits in the case of SPECT.
However, if the reconstruction method considers events one-by-one in their order, then the reconstruction is of a list mode. A list mode reconstruction may involve more information, like the timing of individual events. The differences between binned and list mode reconstruction are given in more detail below.
The required output of a reconstruction method of any type is the emission density function x({right arrow over (v)}) in points {right arrow over (v)} of a volume of interest □, which is typically discretized on a voxel grid. In the case of binned reconstruction the input of the reconstruction method is the measurement data, the measured number of hits in the conceptual detectors, which is connected to the output by a system matrix. An element of the system matrix represents the probability that a particle generated in a voxel is detected by a conceptual detector.
For typical systems, both the number of conceptual detectors and the number of voxels may be in the range of several hundred millions, thus the system matrix has a very large size and the reconstruction method must scale up well and must be appropriate for high performance computation platforms. The known time-consuming processes of iterative tomography reconstruction are targeted by pre-computing system dependent information and always needed more efficient computational platforms (A. J. Reader and H. Zaidi: Advances in PET Image Reconstruction, PET Clinics, Vol. 2, pp. 173-190 (2007)).
Using a pre-computed (U.S. Pat. No. 7,332,721 B2, S. Moehrs, M. Defrise, N. Belcari, A Del Guerra, A. Bartoli, S. Fabbri, and G. Zanetti: Multi-ray-based system matrix generation for 3D PET reconstruction, Phys. Med. Biol., Vol. 53, pp. 6925-6945 (2008) and J. Qi, R. M. Leahy, S. R. Cherry, A. Chatziioannou, and T. H. Farquhar: High-resolution 3D Bayesian image reconstruction using the microPET small-animal scanner. Physics in Medicine and Biology, Vol. 43, pp. 1001-1013 (1998)) or a measured (V. Y. Panin, F. Kehren, C. Michel, and M. Casey: Fully 3-D PET reconstruction with system matrix derived from point source measurements, IEEE Transactions on Medical Imaging, Vol. 25, pp. 907-921 (2006)) system matrix seems to be an attractive option, but in high resolution real three-dimensional scanners this would pose prohibitive memory requirements disadvantageously, even if the matrix is compressed, decomposed as a product of a series of smaller matrices (J. Qi, R. M. Leahy, S. R. Cherry. A. Chatziioannou, and T. H. Farquhar: High-resolution 3D Bayesian image reconstruction using the microPET small-animal scanner, Physics in Medicine and Biology, Vol. 43, pp. 1001-1013 (1998)), or symmetry is exploited (J. L. Herraiz. S. Espaa. S. Garcia, R. Cabido, A. S. Montemayor, M. Desco, J. J. Vaquero, and J. M. Udias: GPU acceleration of a fully 3D iterative reconstruction software for PET using CUDA, in: Nuclear Science Symposium Conference Record (NSS/MIC), 2009 IEEE, pp. 4064-4067 (2009)). On the other hand, in case of a pre-computed or measured system matrix, phenomena like object and material dependent positron range, absorption or scattering cannot be considered, disadvantageously.
Among the high-performance computing possibilities, like Field-Programmable Gate Arrays (FPGA) (M. Leeser, S. Coric, E. Miller, H. Yu, and M. Trepanier: Parallel-beam back projection: an FPGA implementation optimized for medical imaging, in: Proc. Tenth Int. Symposium on FPGA, pp. 217-226 (2002) and US 2008/0095300 A1), multi-CPU systems (D. W. Shattuck, J. Rapela, E. Asma, A. Chatzioannou, J. Qi, and R. M. Leahy: Internet2-based 3D PET image reconstruction using a PC cluster, Physics in Medicine and Biology, Vol. 47, pp. 2785-2795 (2002)), cell processors (M. Kacheiries, M. Knaup, and O. Bockenbach: Hyperfast parallel-beam and conebeam back projection using the CELL general purpose hardware, Medical Physics, Vol. 34. pp. 1474-1486 (2007)), and GPUs (Graphics Processing Units) (F. Xu and K. Mueller: Real-time 3D computed tomographic reconstruction using commodity graphics hardware, Physics in Medicine and Biology, pp. 3405-3419 (2007)), in a comparison to other architectures, the massively parallel GPU has proven to be the most cost-effective platform for tomographic reconstruction tasks (N. Gac, S. Mancini, M. Desvignes, and D. Houzet: High speed 3D tomography on CPU, GPU, and FPGA, EURASIP Journal on Embedded Systems, Article ID 930250 (2008)).
GPUs can be programmed with two different programming models. Shader APIs (Application Programming Interface) present the GPU hardware as the direct implementation of the incremental rendering pipeline (L. Szirmay-Kalos, L. Szcsi. M. Sbert: GPU-Based Techniques for Global Illumination Effects, Morgan and Claypool Publishers, San Rafael, USA (2008)), including both programmable and fixed processing stages. On the other hand, GPGPU (General Purpose GPU) APIs like CUDA (Compute Unified Device Architecture) (NVIDIA, http://developer.nvidia.com/cuda, in: The CUDA Homepage) and OpenCL (Open Computing Language) provide access to the multiprocessors of the GPU where each multiprocessor contains a set of scalar processors sharing the instruction unit, and therefore acting as a SIMD (Single Instruction. Multiple Data) hardware. If only the geometric projection is considered in a tomographic reconstruction, then the shader API implementation is more appropriate since rasterizer and alpha-blending units accessible only through the shader API support these simple calculations (B. Bai, A. Smith: Fast 3D iterative reconstruction of PET images using PC graphics hardware, in: IEEE Nuclear Science Symposium, pp. 2787-2790 (2006); F. Xu and K. Mueller: Real-time 3D computed tomographic reconstruction using commodity graphics hardware, Physics in Medicine and Biology, pp. 3405-3419 (2007); G. Pratx, C. Chinn, P. D. Olcott, C. S. Levin: Fast, accurate and shift-varying line projections for iterative reconstruction using the GPU, IEEE Trans. on Medical Imaging, Vol. 28, pp. 435-445 (2009)). A tomographic imaging technique for GPU is described in F. Xu and K. Mueller: GPU-Acceleration of Attenuation and Scattering Compensation in Emission Computed Tomography, in: 9th Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine '07 (2007). However, when an algorithm gets more complex, the incremental rendering pipeline view of the shader API becomes too restrictive and less intuitive (G. Pratx, L. Xing: GPU computing in medical physics: A review, Med. Phys., Vol. 38, pp. 2685-2698 (2011)), and GPGPU APIs become the better choice for such tasks. The critical issue of the GPGPU programming is thread mapping, i.e. the decomposition of the algorithm to parallel threads that can run efficiently. For example, while simulating particle transport, it is intuitive to mimic how nature works in parallel, and assign parallel computational threads, for example, to randomly generated photon paths (A. Wirth, A. Cserkaszky, B. Kári, D. Légrády, S. Fehér, S. Czifrus, B. Domonkos: Implementation of 3D Monte Carlo PET reconstruction algorithm on GPU, in: Nuclear Science Symposium Conference Record (NSS/MIC), IEEE, pp. 4106-4109 (2009)). However, while significant speed ups can be obtained with respect to a CPU implementation, this somehow natural thread mapping cannot exploit the full potential of GPUs.
In an iterative tomographic reconstruction, forward projection computing the expected detector hits from an actual estimation on emission density data and back projection correcting the emission density data based on the measured and expected detector responses, alternate. Equations of forward projection and back projection are similar in the way that they take many input values—voxel intensities and data of conceptual detectors, respectively—and compute many unknown values—again, data of conceptual detectors and voxel intensities, respectively. This kind of “many to many” computation can be organized in two different ways. Known values can be taken one-by-one obtaining the contribution of a single known value to all of the unknowns, and accumulate the contributions as the different known values are visited. This scheme is called scattering. The opposite—sometimes called orthogonal—approach takes unknown values one-by-one, and obtains the contribution of all known values to this particular unknown value. This approach is called gathering. Generally, input-driven algorithms are of scattering type, while output-driven ones are of gathering type.
A list mode PET reconstruction using GPU architecture is described in US 2011/0182491 A1 and further detailed in a paper (J. Cui, G. Pratx, S. Prevrhal, C. S. Levin: Fully 3D list-mode time-of-flight PET image reconstruction on GPUs using CUDA, Med. Phys., Vol. 38, No. 12 (2011)). These documents state that for a reconstruction using a GPU it is worth if both the forward projection and the back projection is of gathering type. However, the back projection detailed in the documents is of scattering type.
In U.S. Pat. No. 7,778,392 B1 such a computer tomography reconstruction method is described, where the back projection is optimized for GPU architecture. In U.S. Pat. No. 7,876,944 B2 a medical image reconstruction method is described, where a GPU-optimized back projection step is applied.
In U.S. Pat. No. 6,804,325 B1 such a reconstruction method is described, where the end points of the LORs are randomly distributed on the surface of a detector element. In US 2007/0201611 A1 such a PET reconstruction method for a GPU is described, in which tri-linear interpolation is used. In U.S. Pat. No. 7,333,107 B2 a gather-scatter approach is detailed in connection with the forward- and back projection, and the reuse of an already computed data is described for a GPU. In U.S. Pat. No. 7,120,283 B2, US 2010/0266178 A1, U.S. Pat. No. 7,873,236 B2 and US 2011/0051893 A1 medical image processing methods are described for a GPU. In US 2010/0128046 A1 Poisson-disk distribution is used for sampling in a GPU architecture.
Unmatching forward and back projection is described for tomographic reconstruction in G. Zeng, G. Gullberg: Unmatched projector/backprojector pairs in an iterative reconstruction algorithm, IEEE Transactions on Medical Imaging, Vol. 19, pp. 548-555 (2000); R. Guedouar, B. Zarrad: A comparative study between matched and mis-matched projection/back projection pairs used with ASIRT reconstruction method, Nuclear Instruments and Methods in Physics Research, Vol. 619, pp. 225-229 (2010); V.-G. Nguyen, S.-J. Lee, M. N. Lee: GPU-accelerated 3D Bayesian image reconstruction from Compton scattered data, Physics in Medicine and Biology, Vol. 56. pp. 2817 (2011); and N.-Y. Lee, Y. Choi: Theoretical investigation on an unmatched backprojector for iterative reconstruction in emission computed tomography, Journal of the Korean Physical Society, Vol. 59, pp. 367-375 (2011).
In US 2008/0095300 A1, U.S. Pat. Nos. 7,769,217 B2 and 8,000,513 B2 it is described that the forward projection and the back projection are advantageously matching algorithms. It is described in U.S. Pat. No. 7,856,129 B2 that the method detailed therein can be generalized to unmatching projections. In U.S. Pat. Nos. 5,414,623 and 7,596,202 B2 such methods are described, where a filtering is applied between forward and back projection.
In US 2011/0164799 A1 and US 2011/0194747 A1 such methods are described where filtering is applied for the reconstruction data in the Fourier space.
In U.S. Pat. Nos. 7,381,960 B1, 7,759,647 B2 and 7,835,782 B2 the problem of positron propagation before annihilation is detailed. In U.S. Pat. No. 7,835,782 B2 a scheme for handling this problem is described. In U.S. Pat. No. 7,923,690 B2 and US 2011/0248765 A1 it is described that the problem of the positron propagation is not handled.
Tomographic systems are described in U.S. Pat. Nos. 6,373,059 B1, 7,211,799 B2, US 2010/0108896 A1 and US 2011/0127434 A1, in which the different sensitivities of the detector elements on the detection of an event are taken into account. In US 2007/0221855 A1 and WO 2011/036436 A1 the photon transport through adjacent detector elements is handled.
In U.S. Pat. Nos. 7,015,477 B2, 7,332,721 B2, 7,876,941 B2, US 2011/0044546 A1, US 2011/0164031 A1 and U.S. Pat. No. 7,983,465 B2 the possible solutions of the equations of forward- and/or back projection are described.
A tomographic reconstruction method is disclosed in M. Magdics, L. Szirmay-Kalos, Á. Szlavecz, G. Hesz, B. Benyó, Á. Cserkaszky, J. Lantos, D. Légrády, Sz. Czifrus, A. Wirth, B. Kári, G. Patay, D. Vögyes, T. Bükki, P. Major, G. Németh, B. Domonkos: TeraTomo project: a fully 3D GPU based reconstruction code for exploiting the imaging capability of the NanoPET/CT system, 2010 World Molecular Imaging Congress, 11 Sep. 2010. This paper discusses direct Monte Carlo particle transport and presents an adjoint Monte Carlo method which applies gathering type approach both in forward and back projection. Both direct photon hits and scattered photon hits are considered.
Several methods for handling gamma photon scattering in the detector modules are disclosed in M. Magdics, B. Tóth, L. Szécsi, B. Csébfalvi, L. Szirmay-Kalos, Á. Szlavecz, G. Hesz, B. Benyó, Á. Cserkaszky, D. Légrády. Sz. Czifrus, A. Wirth, B. Kári, J. Lantos, G. Patay, D. Völgyes, P. Major. G. Németh, T. Bükki, B. Domonkos: Detector Modelling Techniques for Pre-Clinical 3D PET Reconstruction on the GPU, Proceedings of the 11th International Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine, 11 Jul. 2011, pp. 375-378.
A photon transport simulation method, that executes the very same instruction in parallel threads at a time and eliminates the conditional instructions, is disclosed in L. Szirmay-Kalos, B. Tóth, M. Magdics, D. Légrády, and A. Penzov: Gamma Photon Transport on the GPU for PET, in: Large-scale scientific computing, Springer Berlin, Heidelberg, pp. 435-442, 4 Jun. 2009.
It is a disadvantage of many of the prior art solutions that the level of optimization of the reconstruction methods used for the parallel architectures is not high enough.